Preparations

Load the necessary libraries

library(car)       #for regression diagnostics
library(broom)     #for tidy output
library(broom.mixed)
library(ggfortify) #for model diagnostics
library(sjPlot)    #for outputs
library(knitr)     #for kable
library(effects)   #for partial effects plots
library(ggeffects) #for effects plots in ggplot
library(emmeans)   #for estimating marginal means
library(MASS)      #for glm.nb
library(MuMIn)     #for AICc
library(tidyverse) #for data wrangling
library(DHARMa)   #for residuals and diagnostics
library(nlme)     #for lme
library(lme4)      #for glmer
library(glmmTMB)    #for glmmTMB
library(performance) #for diagnostic plots
library(see)         #for diagnostic plots

Scenario

To investigate synergistic coral defence by mutualist crustaceans, Mckeon et al. (2012) conducted an aquaria experiment in which colonies of a coral species were placed in a tank along with a preditory seastar and one of four symbiont combinations:

  • no symbiont,
  • a crab symbiont
  • a shrimp symbiont
  • both a crab and shrimp symbiont.

The experiments were conducted in a large octagonal flow-through seawater tank that was partitioned into eight sections, which thereby permitted two of each of the four symbiont combinations to be observed concurrently. The tank was left overnight and in the morning, the presence of feeding scars on each coral colony was scored as evidence of predation. The experiments were repeated ten times, each time with fresh coral colonies, seastars and symbiont.

The ten experimental times represent blocks (random effects) within which the symbiont type (fixed effect) are nested.

Read in the data

mckeon = read_csv('../public/data/mckeon.csv', trim_ws=TRUE)
## 
## ── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────
## cols(
##   BLOCK = col_double(),
##   PREDATION = col_double(),
##   SYMBIONT = col_character()
## )
glimpse(mckeon)
## Rows: 80
## Columns: 3
## $ BLOCK     <dbl> 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10…
## $ PREDATION <dbl> 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, …
## $ SYMBIONT  <chr> "none", "none", "none", "none", "none", "none", "none", "non…

Exploratory data analysis

Model formula: \[ y_i \sim{} \mathcal{N}(n, p_i)\\ ln\left(\frac{p_i}{1-p_1}\right) =\boldsymbol{\beta} \bf{X_i} + \boldsymbol{\gamma} \bf{Z_i} \]

where \(\boldsymbol{\beta}\) and \(\boldsymbol{\gamma}\) are vectors of the fixed and random effects parameters respectively and \(\bf{X}\) is the model matrix representing the overall intercept and effects of symbionts on the probability of the colony experiencing predation. \(\bf{Z}\) represents a cell means model matrix for the random intercepts associated with individual coral colonies.

ggplot(mckeon, aes(y=PREDATION, x=SYMBIONT)) +
    geom_point(position=position_jitter(width=0.2, height=0))+
    facet_wrap(~BLOCK)

Fit the model

Model validation

Partial plots

glmer (lme4)

plot_model

plot_model(mckeon.glmer1,  type='eff')
## $SYMBIONT

allEffects

plot(allEffects(mckeon.glmer1))

ggpredict

ggpredict(mckeon.glmer1) %>% plot
## $SYMBIONT

ggemmeans

ggemmeans(mckeon.glmer1,  ~SYMBIONT) %>% plot

glmmTMB (glmmTMB)

plot_model

plot_model(mckeon.glmmTMB1,  type='eff')
## $SYMBIONT

allEffects

plot(allEffects(mckeon.glmmTMB1))

ggpredict

ggpredict(mckeon.glmmTMB1) %>% plot
## $SYMBIONT

ggemmeans

ggemmeans(mckeon.glmmTMB1,  ~SYMBIONT) %>% plot

Model investigation / hypothesis testing

glmer (lme4)

summary

summary(mckeon.glmer1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: PREDATION ~ SYMBIONT + (1 | BLOCK)
##    Data: mckeon
## 
##      AIC      BIC   logLik deviance df.resid 
##     70.7     82.6    -30.4     60.7       75 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -23.2076  -0.1730   0.0976   0.2980   0.8944 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  BLOCK  (Intercept) 11.81    3.437   
## Number of obs: 80, groups:  BLOCK, 10
## 
## Fixed effects:
##                Estimate Std. Error z value Pr(>|z|)   
## (Intercept)       5.096      1.811   2.814  0.00490 **
## SYMBIONTcrabs    -3.842      1.465  -2.623  0.00871 **
## SYMBIONTshrimp   -4.431      1.551  -2.856  0.00429 **
## SYMBIONTboth     -5.599      1.724  -3.247  0.00116 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) SYMBIONTc SYMBIONTs
## SYMBIONTcrb -0.625                    
## SYMBIONTshr -0.652  0.736             
## SYMBIONTbth -0.687  0.732     0.767

Conclusions:

  • the coefficients are presented on a logit scale. Whilst this is not relavant for the purpose of inference testing, it does make it difficult to interpret the coefficients.
  • if we exponentiate the coefficients (\(log(\frac{\rho}{1-\rho})\) -> \(\frac{\rho}{1-\rho}\)), they will be presented on a odds ratio scale, and thus:
    • the intercept (none symbionts) will be 163.33. That is, corals without a symbiont are 163.33 times more likely to be preditated on than not predated on. The odds of predation in this the absence of symbionts is 163.33:1.
    • in the presence of a crab symbiont, the odds of being predated on are only 0.02 times that of the none symbiont group. That is, in the presence of a crab symbiont, the odds of predation decline by 98%.
    • in the presence of a shrimp symbiont, the odds of being predated on are only 0.01 times that of the none symbiont group. That is, in the presence of a shrimp symbiont, the odds of predation decline by 99%.
    • in the presence of both crab and shrimp symbionts, the odds of being predated on are only 0 times that of the none symbiont group. That is, in the presence of both crab and shrimp symbiont, the odds of predation decline by 100%.
  • if we backtransform the intercept full to the response scale (probability scale), (\(log(\frac{\rho}{1-\rho})\) -> \(\rho\)), the intercept is interpreted as the probability that corals will be predated in the absence of of symbionts is 1

tidy

tidy(mckeon.glmer1, effect='fixed', conf.int=TRUE, infer=c(TRUE, TRUE))
tidy(mckeon.glmer1, effect='fixed', conf.int=TRUE, infer=c(TRUE, TRUE),  exponentiate=TRUE)
tidy(mckeon.glmer1, effect='fixed', conf.int=TRUE, infer=c(TRUE, TRUE),  exponentiate=TRUE) %>% kable
effect term estimate std.error statistic p.value conf.low conf.high
fixed (Intercept) 163.3256181 295.7986843 2.813623 0.0048987 4.6929398 5684.1252579
fixed SYMBIONTcrabs 0.0214459 0.0314110 -2.623280 0.0087088 0.0012151 0.3785015
fixed SYMBIONTshrimp 0.0119023 0.0184657 -2.856063 0.0042893 0.0005689 0.2490134
fixed SYMBIONTboth 0.0037016 0.0063822 -3.247315 0.0011650 0.0001261 0.1086478

Conclusions:

  • the coefficients are presented on a logit scale. Whilst this is not relavant for the purpose of inference testing, it does make it difficult to interpret the coefficients.
  • if we exponentiate the coefficients (\(log(\frac{\rho}{1-\rho})\) -> \(\frac{\rho}{1-\rho}\)), they will be presented on a odds ratio scale, and thus:
    • the intercept (none symbionts) will be 163.33. That is, corals without a symbiont are 163.33 times more likely to be preditated on than not predated on. The odds of predation in this the absence of symbionts is 163.33:1.
    • in the presence of a crab symbiont, the odds of being predated on are only 0.02 times that of the none symbiont group. That is, in the presence of a crab symbiont, the odds of predation decline by 98%.
    • in the presence of a shrimp symbiont, the odds of being predated on are only 0.01 times that of the none symbiont group. That is, in the presence of a shrimp symbiont, the odds of predation decline by 99%.
    • in the presence of both crab and shrimp symbionts, the odds of being predated on are only 0 times that of the none symbiont group. That is, in the presence of both crab and shrimp symbiont, the odds of predation decline by 100%.
  • if we backtransform the intercept full to the response scale (probability scale), (\(log(\frac{\rho}{1-\rho})\) -> \(\rho\)), the intercept is interpreted as the probability that corals will be predated in the absence of of symbionts is 1

tab_model

# warning this is only appropriate for html output
sjPlot::tab_model(mckeon.glmer1, show.se=TRUE, show.aic=TRUE)
  PREDATION
Predictors Odds Ratios std. Error CI p
(Intercept) 163.33 295.80 4.69 – 5684.13 0.005
SYMBIONT [crabs] 0.02 0.03 0.00 – 0.38 0.009
SYMBIONT [shrimp] 0.01 0.02 0.00 – 0.25 0.004
SYMBIONT [both] 0.00 0.01 0.00 – 0.11 0.001
Random Effects
σ2 3.29
τ00 BLOCK 11.81
ICC 0.78
N BLOCK 10
Observations 80
Marginal R2 / Conditional R2 0.228 / 0.832
AIC 70.706

glmmTMB (glmmTMB)

summary

summary(mckeon.glmmTMB1)
##  Family: binomial  ( logit )
## Formula:          PREDATION ~ SYMBIONT + (1 | BLOCK)
## Data: mckeon
## 
##      AIC      BIC   logLik deviance df.resid 
##     62.9     74.8    -26.4     52.9       79 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  BLOCK  (Intercept) 17.83    4.223   
## Number of obs: 80, groups:  BLOCK, 10
## 
## Conditional model:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       4.420      1.782   2.481 0.013109 *  
## SYMBIONTcrabs    -3.317      1.313  -2.526 0.011534 *  
## SYMBIONTshrimp   -3.956      1.416  -2.793 0.005218 ** 
## SYMBIONTboth     -5.152      1.565  -3.291 0.000997 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Conclusions:

  • the coefficients are presented on a logit scale. Whilst this is not relavant for the purpose of inference testing, it does make it difficult to interpret the coefficients.
  • if we exponentiate the coefficients (\(log(\frac{\rho}{1-\rho})\) -> \(\frac{\rho}{1-\rho}\)), they will be presented on a odds ratio scale, and thus:
    • the intercept (none symbionts) will be 83.1. That is, corals without a symbiont are 83.1 times more likely to be preditated on than not predated on. The odds of predation in this the absence of symbionts is 83.1:1.
    • in the presence of a crab symbiont, the odds of being predated on are only 0.04 times that of the none symbiont group. That is, in the presence of a crab symbiont, the odds of predation decline by 96%.
    • in the presence of a shrimp symbiont, the odds of being predated on are only 0.02 times that of the none symbiont group. That is, in the presence of a shrimp symbiont, the odds of predation decline by 98%.
    • in the presence of both crab and shrimp symbionts, the odds of being predated on are only 0.01 times that of the none symbiont group. That is, in the presence of both crab and shrimp symbiont, the odds of predation decline by 99%.
  • if we backtransform the intercept full to the response scale (probability scale), (\(log(\frac{\rho}{1-\rho})\) -> \(\rho\)), the intercept is interpreted as the probability that corals will be predated in the absence of of symbionts is 0.99

tidy

tidy(mckeon.glmmTMB1, effect='fixed', conf.int=TRUE, infer=c(TRUE, TRUE))
## Warning in confint.glmmTMB(x, method = tolower(conf.method), level =
## conf.level, : extra arguments ignored: infer
tidy(mckeon.glmmTMB1, effect='fixed', conf.int=TRUE, infer=c(TRUE, TRUE),  exponentiate=TRUE)
## Warning in confint.glmmTMB(x, method = tolower(conf.method), level =
## conf.level, : extra arguments ignored: infer

Conclusions:

  • the coefficients are presented on a logit scale. Whilst this is not relavant for the purpose of inference testing, it does make it difficult to interpret the coefficients.
  • if we exponentiate the coefficients (\(log(\frac{\rho}{1-\rho})\) -> \(\frac{\rho}{1-\rho}\)), they will be presented on a odds ratio scale, and thus:
    • the intercept (none symbionts) will be 83.1. That is, corals without a symbiont are 83.1 times more likely to be preditated on than not predated on. The odds of predation in this the absence of symbionts is 83.1:1.
    • in the presence of a crab symbiont, the odds of being predated on are only 0.04 times that of the none symbiont group. That is, in the presence of a crab symbiont, the odds of predation decline by 96%.
    • in the presence of a shrimp symbiont, the odds of being predated on are only 0.02 times that of the none symbiont group. That is, in the presence of a shrimp symbiont, the odds of predation decline by 98%.
    • in the presence of both crab and shrimp symbionts, the odds of being predated on are only 0.01 times that of the none symbiont group. That is, in the presence of both crab and shrimp symbiont, the odds of predation decline by 99%.
  • if we backtransform the intercept full to the response scale (probability scale), (\(log(\frac{\rho}{1-\rho})\) -> \(\rho\)), the intercept is interpreted as the probability that corals will be predated in the absence of of symbionts is 1

tab_model

# warning this is only appropriate for html output
sjPlot::tab_model(mckeon.glmmTMB1, show.se=TRUE, show.aic=TRUE)
  PREDATION
Predictors Odds Ratios std. Error CI p
(Intercept) 83.10 148.05 2.53 – 2730.03 0.013
SYMBIONT [crabs] 0.04 0.05 0.00 – 0.48 0.012
SYMBIONT [shrimp] 0.02 0.03 0.00 – 0.31 0.005
SYMBIONT [both] 0.01 0.01 0.00 – 0.12 0.001
Random Effects
σ2 3.29
τ00 BLOCK 17.83
ICC 0.84
N BLOCK 10
Observations 80
Marginal R2 / Conditional R2 0.149 / 0.867
AIC 62.861

Further analyses

Summary figure

glmer (lme4)

emmeans(mckeon.glmer1, ~SYMBIONT, type='response') %>%
  as.data.frame %>%
  ggplot(aes(y=prob,  x=SYMBIONT)) +
  geom_pointrange(aes(ymin=asymp.LCL,  ymax=asymp.UCL))

glmmTMB (glmmTMB)

emmeans(mckeon.glmmTMB1, ~SYMBIONT, type='response') %>% 
  as.data.frame %>%
  ggplot(aes(y=prob,  x=SYMBIONT)) +
  geom_pointrange(aes(ymin=lower.CL,  ymax=upper.CL))

# References

Mckeon, Seabird, Adrian Stier, Shelby Mcilroy, and Benjamin Bolker. 2012. “Multiple Defender Effects: Synergistic Coral Defense by Mutualist Crustaceans.” Oecologia 169 (February): 1095–1103. https://doi.org/10.1007/s00442-012-2275-2.